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We introduce the concept of concatenated tensor networks to efficiently describe quantum states. 
We show that the corresponding concatenated tensor network states can efficiently describe time 
evolution and possess arbitrary block-wise entanglement and long-ranged correlations. We illustrate 
the approach for the enhancement of matrix product states, i.e. ID tensor networks, where we 
replace each of the matrices of the original matrix product state with another ID tensor network. 
This procedure yields a 2D tensor network, which includes - already for tensor dimension two - 
all states that can be prepared by circuits of polynomially many (possibly non-unitary) two-qubit 
quantum operations, as well as states resulting from time evolution with respect to Hamiltonians 
with short-ranged interactions. We investigate the possibility to efficiently extract information from 
these states, which serves as the basic step in a variational optimization procedure. To this aim 
we utilize known exact and approximate methods for 2D tensor networks and demonstrate some 
improvements thereof, which are also applicable e.g. in the context of 2D projected entangled pair 
states. We generalize the approach to higher dimensional- and tree tensor networks. 

PACS numbers: 03.67.Mn,03.65.Ud,03.67.Lx,02.70.-c 



Contents 



B. Monte-Carlo sampling of CTS 



I. Introduction 

II. Concatenated tensor networks 

A. MPS and PEPS 

B. Concatenated tensor network states 



1 

2 
2 
3 



III. Properties of concatenated tensor network 
states 

A. Concatenated MPS 

1. Interpretation in terms of (post selected) 
quantum circuits 

2. Description of time evolution 

3. Interpretation in terms of measurement-based 
quantum computation. 

4. Interpretation as MPS with large effective 
dimension 

B. Concatenated PEPS 

C. Tensor tree states with internal structure 

IV. Applications 

A. The descriptive potential of CTS 

B. Reading out information from CTS 

1. Approximate contraction of 2D tensor 
networks 

2. An error correction scheme 

3. Compressibility of sequences of matrix 
product operators 

4. Special cases - exact contraction 



10 
11 

5. The advantage of CTS in the operator picturdll 
Summary and outlook [ill 
Acknowledgments [l2l 



References 



I. INTRODUCTION 

The classical simulation of complex quantum systems is 
one of the central problems in modern physics. Given the 
exponential growth of the state space with the system size, 
such a classical simulation seems infeasible. However, it has 
been realized that quantum systems occurring in nature of- 
ten do only populate a small subspace. Identifying this sub- 
space is hence the first step towards a successful classical 
simulation. For ground states of (non-critical) strongly cor- 
related quantum spins in a one dimensional setup, matrix- 
product states (MPS) 0, 0, 0] turn out to provide a proper 
parametrization for this subspace 0, MPS can not 

only efficiently describe such ground states, but it is also 
possible to efficiently read out physical information from 
this description, e.g., to compute expectation values of lo- 
cal observables and correlation functions. Moreover, MPS 
form the basis of the density matrix renormalization group 
(DMRG) [HQ; a powerful numerical method that has been 
successfully applied to various problems in ID. The relation 
between the DMRG and MPS is an example how physi- 
cal insight into the logic of a preparation (renormalization) 
procedure can be manifestly encoded into the structure of 
a state class. 

Recent approaches to simulate ground states of strongly 
correlated systems in critical systems or higher dimensions 
follow a similar approach. A variety of states such as pro- 
jected entangled pairs (PEPS) |2fll . sequentially generated 
stat es jla |. string-bond states [33], weighted graph states 
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23, 24], renormalization ansatz with graph enhancement 
251 ] or the multiscale entanglement renormalization ansatz 
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(MERA) [Hi have been introduced with the aim of effi- 
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ciently parametrizing the relevant subspace. The entan- 
glement properties of the corresponding states form the 
guideline and determine the potential applicability of the 
methods. For instance, MERA can provide a logarithmic 
divergence for block-wise entanglement in critical ID sys- 
tems, while e.g. 2D variants of MERA as well as PEPS and 
string-bond states fulfill area laws for block-wise entangle- 
ment, typically to be found in ground states of 2D systems. 
In all cases it is crucial that not only an efficient description 
of the states can be obtained, but also that information can 
be efficiently extracted, either exactly or in an approximate 
way. Based on these states, variational methods for ground 
state approximation and (real and imaginary) time evolu- 
tion have been developed and tested. While MPS, MERA 
and PEPS lead to good descriptions of ground states for 
non-critical ID systems, critical ID systems and 2D sys- 
tems respectively, none of the proposed classes seems to 
be suitable to properly describe time evolution. In fact, it 
has been argued that simulating time evolution is in gen- 
eral hard [43|, as the block- wise entanglement grows - al- 
ready for ID systems - linearly in time, leading quickly to 
a volume law. The entanglement contained in an MPS is 
bounded by the dimension of the matrices or tensors, and 
the entanglement contained in a PEPS follows an area law. 

Here we present a class of tensor network states for which 
such limitations do not apply, and which allow one in prin- 
ciple to efficiently describe states resulting from time evo- 
lution or quantum computation. To construct these states 
we make use of the basic idea underlying previous tensor 
network structures. In these structures, a simplification of 
the existing description can be achieved by replacing ten- 
sors of high rank (i.e., with many indices) by a network of 
tensors of low rank (i.e., with few indices) with appropriate 
topology. The choice of the underlying tensor network de- 
termines qualitatively different sub-classes of states, in pre- 
vious approaches e.g. having lead to MPS or PEPS when 
describing a ID or 2D structure respectively. We apply this 
idea in an iterative, or concatenated, fashion, leading to con- 
catenated tensor network states (CTS). That is, each of the 
tensors appearing in a tensor network is itself repeatedly re- 
placed by another tensor network. The resulting structure 
is again a tensor network, similar to a PEPS, with the main 
difference that only some of the tensors are associated with 
physical particles. 

The efficient and exact extraction of information, e.g., ex- 
pectation values or correlation functions, from an arbitrary 
tensor network is in general not possible, as they rely on a 
contraction of the network, i.e., summations over all indices 
of the network. Even for 2D tensor networks, the contrac- 
tion is known to be computationally hard (#P-hard) [42]. 
However, for certain special cases exact evaluation is possi- 
ble. In addition, also approximate contraction and certain 
Monte-Carlo methods have been developed and successfully 
applied in the context of 2D PEPS and imaginary time 
evolution [29]. We demonstrate the applicability of the es- 
tablished methods to the CTS and several enhancements 
thereof. We moreover demonstrate that there are novel im- 
plementations of algorithms like (imaginary) time evolution 
of ID systems and the application of quantum circuits that 



are more efficient in the CTS than in MPS. 

This paper is structured as follows. In Sec. HU we will 
introduce the CTS, give examples and illustrate their prop- 
erties from an analytic point of view. In Sec. IIIIl we discuss 
the applications of CTS and illustrate the potential of the 
CTS to describe states relevant in physics. As an exam- 
ple, we give the numerical treatment of a toy model, more 
precisely, we will describe a state originating from the time 
evolution of a product state governed by the Ising Hamil- 
tonian. In Sec. IIVI we finally show several ways to extract 
information from a CTS, thereby utilizing and improving 
methods to (approximately) contract 2D tensor networks. 

II. CONCATENATED TENSOR NETWORKS 

In this section we introduce the CTS and in the context 
of the problems having lead to tensor network descriptions 
in general. 

A generic quantum state of N d-level systems can be 
written in a basis whose elements are tensor products of 
basis states of the local d-level systems. The quantum state 
is then characterized by the coefficients of these basis states, 
which are tensors A SlS2 ,,, SN of rank N and dimension d 

d 

\ip)= A SlS2 ... SN \ Sl s 2 ■ ■ -s N ). (1) 

Sl,S2,...,S n =l 

Hence, the description of such a state consists of d N com- 
plex parameters. This exponential growth of the number 
of parameters used in the generic description makes it un- 
suitable for numerical analysis. 

A. MPS and PEPS 

The tensor A SlS2 ... SN of the generic description given 
above can be decomposed into a tensor network, thereby 
imposing a structure in this set of parameters. To do so, we 
will represent A S1S2 ... SN by another set of tensors of smaller 
rank. Some of the indices of the small-rank tensors corre- 
spond to the state of a physical site {si, s 2 ■ ■ ■ ,sn} as be- 
fore. The remaining auxiliary indices are shared between 
pairs of the small-rank tensors, and to recover the coeffi- 
cient of a basis state of the physical system, the shared 
indices will be contracted, i.e., summed over. The infor- 
mation which tensors share indices can be represented by a 
graph, where "tensors" correspond to vertices and "sharing 
an index" corresponds to an edge. The indices correspond- 
ing to physical states will in the following be called open. 
We will furthermore use greek letters ctj, etc. to refer to 
shared indices, while open indices will be denoted by sj. 

As an example, one may use a ID structure for the tensor 
network, leading to MPS (see Fig. QJi) 

A 

D 

= A[H] A[*2] A[SN-l] A[sn] (r)\ 

/ j cti ^-Q2Ct3 • ' ctiv-iajv Qw ' \ ) 
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FIG. 1: (Color online), (a) Graphical representation of a ID 
tensor network (MPS). The boxes correspond to tensors, where 
shared indices are summed over. Open indices correspond to 
physical particles (red tensors), (b) Each of the tensors in the 
original tensor network is replaced by a ID tensor network (ma- 
trix product operator) arranged in y-direction. Auxiliary ten- 
sors (no open indices) are drawn in blue. This leads to a 2D 
tensor network. 

which are described by the tensors Aa*a i+1 and ^4a i ^ ) ^a^v ^ ■ 
For a fixed choice of S1S2 . . . sn, the coefficient A SlS2 ,,, SN is 
obtained by calculating the product of the D x D matrices 
Aq.'L-h (except at the border, where one has vectors Aa^). 
By choosing D large enough (but still D < d N ), one can 
represent any tensor and hence any quantum state in this 
form. A restriction to small D allows to describe a certain 
subset of states efficiently. 

In a similar way, one can consider tensor network struc- 
tures with different topology and higher dimensional con- 
nectivity. If the physical system consists of particles on 
a 2D regular lattice, the procedure analogous to the con- 
struction of the MPS described above yields a 2D regular 
tensor grid, e.g., 

D 

A A^ A [S21 A^ A^ 

^s 1S2 ...s 9 — /1 Qi,ft /1 aia 2 A a 2 /3 3 ^/3ia 3 /34 

greek indices— 1 

v A^ A^®} A^ S7 ^ A^ s ^ 4^°] 

the structure corresponding to the projected entangled pair 
states (PEPS) in 2D (see Fig. EK). 

Tensor networks that have been subject to detailed in- 
vestigation include one-dimensional graphs with and with- 
out periodic boundary conditions (MPS) , trees (3ll [33| and 
two-dimensional lattices (PEPS) [2g]. Investigations of net- 
works of different topology have shown that ID and tree- 
like structures are generally easy to simulate numerically. 
Tensor networks corresponding to graphs with man y lo ops, 
on the other hand, are generally hard to simulate [3ll . l33| 
and only in special cases efficient algorithms are known, 
see, e.g., Ref. [21(. Some of the networks, e.g. correspond- 
ing to a 2D lattice, are even known to correspond to states 
being resources of measurement based quantum computa- 
tion and hence (having a generally applicable method) to 
treat these tensor network states numerically efficiently and 
exactly would mean to efficiently simulate a quantum com- 
puter classically. In fact, the contraction of such 2D net- 
works was proven to be a computationally hard problem in 
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FIG. 2: (Color online). Examples of concatenated 2D tensor 
network states. The boxes correspond to tensors, where joint 
indices are summed over. Open indices correspond to physical 
particles (red tensors), while auxiliary tensors (no open indices) 
are drawn in blue, (a) Original 2D tensor network, where each 
of the tensors corresponds to a physical particle, (b) Each of 
the original tensors is replaced by a ID tensor network (MPS, 
consisting of 3 tensors, two of which are auxiliary tensors) in 
horizontal direction, (c) Each of the original tensors is replaced 
by a 2D tensor network (of size 3x3) arranged in the same plane 
as the original 2D tensor network, (d) Each of the initial tensors 
is replaced by an MPS perpendicular to the original plane (z- 
direction). This leads a 3D tensor network structure. 

general [52j. 



B. Concatenated tensor network states 

We will consider concatenated tensor networks in the fol- 
lowing. That is, given a tensor network as in the previous 
subsections, we will replace each individual tensor in the 
network by another tensor network. This can in principle 
be done in an iterative way, leading to concatenated tensor 
structures. We will typically only consider tensor networks 
stemming from few iterations, given the fact that the total 
number of tensors increases exponentially with the number 
of iterations. Notice that most of the additional tensors 
that we introduce will be auxiliary tensors, i.e., without 
open indices and hence not corresponding to quantum sys- 
tems. We also remark that it is not necessary to use the 
same tensor structure at each concatenation level. 

The structure that we finally obtain is again a (possi- 
bly high-dimensional) tensor network. As long as the total 
number of tensors, as well as their rank and dimension, is 
polynomially bounded, we obtain a class of states that can 
be described by a polynomial number of parameters, i.e., 
efficiently. We call the family of quantum states that can be 
described in this way concatenated tensor networks states 
(CTS). 

The key element of this approach is to impose internal 
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structure on the tensor description being used, thereby re- 
ducing the information content, while its ability to describe 
entanglement is in principle kept. This allows one to de- 
scribe states with a large amount of block-wise entangle- 
ment, up to a volume law, and long-ranged correlations 
using only small rank tensors of small dimension at the 
elementary level. 

By construction, the imposed structure is similar to the 
one behind the very successful DMRG renormalization. In- 
dependently of this ancestry of the ansatz, there are some 
illuminating interpretations going beyond the DMRG pic- 
ture. Different from the DMRG, the renormalization struc- 
ture in the concatenated tensor network is not necessarily 
applied in a spacial fashion, but (being subject to interpre- 
tation and depending on the actual network) in a timely 
fashion, e.g., as preparatory applications of certain oper- 
ators in a Suzuki- Trotter expansion, or going further, as 
state-preparing applications of generic operators. We find 
that already with a 2D network of size poly(TV) and ten- 
sor dimension two, all states that can be prepared by a 
polynomially sized quantum circuit can be represented as 
CTS. Furthermore, the picture that the CTS stems from a 
preparation using measurement based quantum computa- 
tion (MQC) is possible. All these interpretations are suited 
to inspire further development and nurture some hope that 
the described state class might - by virtue of its construc- 
tion - be suited for a good description of time evolved states 
or quantum circuits. 



III. PROPERTIES OF CONCATENATED TENSOR 
NETWORK STATES 



individual tensors by ID tensor networks can be iterated. 
At the next level, one obtains a 3D tensor network and 
so forth. We remark that one may also consider 2D ten- 
sor networks with periodic boundary conditions, either in 
horizontal or vertical direction. 

In the following we will consider a 2D tensor networks 
(i.e., only the first iteration level) of size N x M with M = 
poly (TV). We analyze the states that can be described by 
such a CTS, and study their entanglement features. We 
show that 

• All states that can be created by a polynomially sized 
quantum circuit can be efficiently described by such a 
2D CTS with Dt = 2. This includes unitary quantum 
circuits as well as post-selected quantum circuits. 

• All states resulting from a time evolution for a time 
t with respect to short-range Hamiltonians can be 
efficiently described by an N x M 2D CTS, where M 
scales quadratically with time t. 



• A subclass of matrix product states with an effective 
bond-dimension of the order of x 
scribed efficiently by an N x M 2D CTS. 



can be de- 



Regarding the entanglement features, we show 

• The block-wise entanglement of an N x M 2D CTS 
can be 0{M). In particular, states with a volume 
law for block- wise entanglement and with long-ranged 
correlations can be described efficiently. 



Interpretation in terms of (post selected) quantum circuits 



In the following, we give a number of examples of CTS 
and discuss their properties. 



A. Concatenated MPS 

We now consider concatenated MPS. We start with a 
ID tensor network as shown in Fig. QJi, and replace each 
of the tensors Aa k a k+1 by a ID tensor network, as shown 

in Fig. [TJd. More precisely, each matrix Aa£a k+1 for Sk — 
1, 2 . . . , d is replaced by a matrix product operator (MPO) 

H, 



A [sk] <-> 

oi k a k+1 

E 



(3) 



.4 



8*] 



and the indices otk are replaced by oP k € (1, 2, . . . , D) corre- 
sponding to several connections to the neighboring tensors. 
Notice that the effective dimension of all these connections 
together is given by x = Yik ^ n tms wa y we obtain a 
2D tensor network, where only N tensors A^^ have open 
indices and correspond to physical sites, while there are 
(N — 1)M auxiliary tensors B. The process of replacing 



Here we show that for a specific choice of tensors the 2D 
tensor network can be interpreted as a quantum circuit con- 
sisting of generic gates. We consider a quantum circuit for 
N qubits of depth M = 0(poly(A r )). We find that one can 
describe the resulting state from such a quantum computa- 
tion by a 2D tensor network of size of order 0(N x M), i.e., 
of polynomially many tensors, where the tensor dimension 
is D = 2, Let us now demonstrate how a standard quan- 
tum circuit consisting of arbitrary single-qubit rotations 
and two-qubit phase gates - which constitute a universal 

gate set - can be encoded into the tensor network. We de- 
fa i) 

note the auxiliary tensors by Bala r a u ad an d the ones con- 
nected to physical particles by Aaf2ra u s d (typically), where 
the sub-indices l,r,u,d stand for left, right, up and down, 
and i,j are labels that indicate the position of the tensor 
in the 2D tensor network (i th row and j th column). The 

uppermost line of tensors Ba l a r0ld have no "up" index, and 
similarly the tensors at the border do not have left/right 
indices. We identify each horizontal line of tensors with a 
certain time step in the circuit, and the first (uppermost) 
line is used to initialize the input state to \Q)® N (or some 
other product state), while the last line corresponds to the 
output state. 

Initialization can, e.g., be achieved by choosing -Bqoo = 
1 and all other entries 0, where we identify the component 
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(1) of the down link with the state |0) (|1)). The basic idea 
is then to either erase the left-right links between two neigh- 
boring tensors, so that processing of individual qubits can 
be performed, or make use of this link to perform an (entan- 
gling) two qubit gate. In the contraction of the tensor net- 
work, one sums over all possible values for each of the links. 
Hence if we choose Va r a u ad : Bq 1 ^ = 0, the link to the 
left is essentially broken [H3| • Similarly, the link to the right 
can be broken by choosing \/aia u ad ■ B ai o auCtd (i, j) = 0. 
Hence the choice 



B 



lla u a d 



(4) 



(and all other entries are 0) allows us to implement the 
single-qubit (unitary) operation 



U= ^2 U aaa Ja d ) (a u \ 



(5) 



on qubit j in time step i, 

For a two-qubit phase gate diagQl, 1,1,-1]), i.e., 

l 

£/pg = Ua 

d t3 d a lt pja! d f3 d ){a u p u \ 

1 

= E (-V ad ' 0d \<Xd0d){a d f3 d \, (6) 

a d ,0 d =O 



acting on qubits j, j + 1 in time step i, we find that the fol- 
lowing choice of tensors allows one implement this gate: 

= 1; 



^1000 — ^lOll 



- #1111 



0100 — -"0111 — 

l(^ini" = — 2, while all other tensors are zero. This 
can be seen by noting that the links to left (particle j — 1) 
and right (particle j + 1) are broken, and by contracting 
the two tensors over their joint index (a r ,f3i). Other two 
qubit gates corresponding to the class of CNOT and phase 
gates [45] (i.e., gates that can create only Schmidt-rank two 
states or only one e-bit of entanglement) can be realized. 
Among these gates are e.g. controlled phase gates with a 
controllable phase <p, Upg(<p) = diag([l, 1, 1, e lip ]). 

To give an example for a subclass of states with a large 
amount of entanglement to be created by operators and to 
be hold by a simple CTS description, consider controlled 
phase gates Upcif) between arbitrary pairs of particles ini- 
tially prepared in |+) = l/\/2(|0) + |1)V These circuits 
prepare weighted graph states (WGS) [13, US], utilizing 
only 0{N 2 ) gates. Using nearest neighbor gates, one needs 
at most 0(N 3 ) phase gates to prepare an arbitrary WGS, 
although one is not restricted to these in our setup. As 
demonstrated in [23], WGS can have maximal block- wise 
entanglement, maximal localizable entanglement as well as 
long-ranged correlations. Similarly, as shown in [4^, typi- 
cal states with O(L) block-wise entanglement for all blocks 
of length L can be generated by 0(N 3 ) two-qubit gates 
acting on arbitrary pairs of particles, leading to a tensor 
network of size A x 0(N 4 ). 

The generalization to other (non-unitary) circuits or 
other elementary gates is straightforward. For instance, the 



unitary matrix U ad0lu in Eq. [4] can be replaced by an arbi- 
trary matrix A adau , corresponding to an arbitrary single- 
qubit operation. In particular, a single-qubit measurement 
with a selected outcome can be described in this way by 
choosing A to be a ID projector. Using such a construc- 
tion, one obtains all states that can be described by an 
arbitrary post-selected quantum circuit. The correspond- 
ing complexity class is postBQP, which is in fact equivalent 
to PP 0]. 

Finally, we remark that, when considering a 2D tensor 
network on a tilted lattice, one can interpret the tensors 
directly as (unitary or non-unitary) quantum gates acting 
on nearest neighbors (see also Ref. 



2. Description of time evolution 

Similarly to the description of a polynomially sized quan- 
tum circuit, one can find, as a special case, a description 
of time evolution in terms of a polynomially sized 2D ten- 
sor network. Consider for example a nearest-neighbor ID 
Hamiltonian H = • -Hjj+i that we decompose into two 
parts, Hi and H 2 , where Hi [H 2 ] contains pairwise com- 
muting terms acting on different systems. That is, Hi = 
J2k H 2k-i,2k, while H 2 = Y,k H M,2k+i, see Refs. @, 0]. 
Using the Suzuki- Trotter expansion, we can write 



e ~itH _ ( ,i(H 1 +H 2 )t 



lim T]{e- lHl 



t/n —iH 2 t/r 



fe=l 



where for a fixed time t we obtain a proper approximation 
with bounded error e by choosing n — 0(t 2 /e), see Ref. [Hll ]. 
and hence a fixed small time step St = t/n = 0(e/t), Hence 
the time evolution for time t is accurately described by a 
sequence of 2n gates of the form e~ iStHj , where n scales 
quadratically with t [13] • Each of the gates e~ StHj , j = 1, 2 
can be described by a 2D tensor network of size Axe, 
where c is a small constant, similarly as discussed for poly- 
nomially sized quantum circuits in the previous subsection. 
The state resulting from a time evolution for time t with 
respect to the Hamiltonian H applied to some initial prod- 
uct state can hence be described by a 2D tensor network of 
size N x M with M = 2cn = 0(t 2 /e). 



3. Interpretation in terms of measurement-based quantum 
computation. 

Another interpretation of such a tensor network descrip- 
tion is provided by measurement-based quantum compu- 
tation (MQC) [H, [H]. One can view the 2D tensor net- 
work as the PEPS description of e.g. a 2D cluster state, 
where all but A^ particles (last row) are measured out. The 
choice of tensors allows one to choose the measurement di- 
rections of the corresponding (auxiliary) particles. In turn, 
the measurement pattern (i.e., the choice of measurements) 
determines the quantum state that is generated at the out- 
put qubits (corresponding to the open legs in our tensor 
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network). In fact, as each choice of tensor corresponds to 
a specific measurement outcome, we consider only a single 
branch of the measurement-based quantum computation, 
i.e. probabilistic MQC with some non-zero success prob- 
ability [i^l . Again, this is equivalent to all post-selected 
quantum circuits. Notice that also other tensor structures 
are universal in this probabilistic sense i.e., allow one 
to describe/generate all quantum states. 

In other words, the tensor network describes a quantum 
state of N + M particles, where the M auxiliary particles 
are measured out in order to finally generate a state of 
N quantum particles. The auxiliary particles (auxiliary 
tensors) allow one to assist the generation of an enlarged 
class of states. 



4- Interpretation as MPS with large effective dimension 

A general MPS corresponding to a ID tensor network 
with matrix dimension x ls described by NO(x 2 ) parame- 
ters. The block-wise entanglement in such a MPS is limited 
by log 2 X- For a 2D CTS of size N x M , and tensors of 
dimension D, we observe that one may still interpret the 
resulting state as an MPS or ID network (by contracting 
the MPO along the vertical direction). The effective ma- 
trix dimension of the corresponding MPS is now given by 
X = D M . This also implies that the potential block- wise 
entanglement, measured by the entropy, between systems 
(1 ... k) and (k - 1 ... TV) is given by log 2 D M = M log 2 D. 
This corresponds to an exponential increase in effective 
bond-dimension while increasing the total number of pa- 
rameters to describe the state only polynomially. Clearly, 
only a specific subset of states with a given block-wise en- 
tanglement can be described by such a 2D tensor network, 
however this set now includes states with large block-wise 
entanglement. If M — O(N), it follows that the corre- 
sponding states can even be maximally entangled, i.e., ful- 
fill a volume law. 

Notice that describing states in terms of such a 2D CTS 
can already be useful for small M. Consider for instance 
ground states of ID critical systems, where it is known that 
a good description in terms of an MPS requires a matrix 
dimension x = 0(2 logN ) [|, Eo[. Similarly, the states re- 
sulting from a time evolution for a time t with respect to a 
nearest-neighbor Hamiltonian possess block-wise entangle- 
ment growing linearly with t, leading eventually to volume 
laws. This implies that a description in terms of a general 
MPS requires matrices of dimension x — 0(2 ), i.e., expo- 
nentially many parameters. In turn, the 2D CTS can pos- 
sess block- wise entanglement scaling as O(MlogD), while 
the total number of parameters is of order 0(MND A ). 
That is, already for D fixed and M = 0(N) a volume law 
can be obtained. For a specific example for the successful 
application of such a CTS description in the context of time 
evolution, see Sec. IIV1 

The natural limitation of the entanglement we describe is 
not given by its quantity, which can be expressed, for exam- 
ple, as the cardinality of the set of Schmidt coefficients in 
a bipartition of the given state. The limitation underlying 



the efficiency is rather introduced by a certain structure, 
or order, within this (potentially very large) set of Schmidt 
coefficients. Depending on the situation, the (itself vari- 
able) structure of the entanglement will not have such a 
big impact on accuracy that the limitation of the quantity 
would have. 



B. Concatenated PEPS 

We now turn to (the CTS extension of) 2D tensor net- 
works of size N x N , or equivalently 2D PEPS. In contrast 
to 2D networks considered in the previous subsection, all 
tensors in such a 2D tensor network have open indices and 
are hence associated with a physical system. As before, 

Is- ■] 

we now replace each of these tensors A a ^a r a u a d by another 
tensor network. There are several possibilities to do this 
(see Fig.[2]), 

(i) We use a ID tensor network (matrix-product oper- 
ator) of dimension D with M tensors, arranged in 
horizontal direction. One of the tensors has an open 
index corresponding to a physical system, while M — 1 
are auxiliary tensors. This leads to a {NM) x N 2D 
tensor network depicted in Fig. (2}d. Similar, one can 
use a ID network arranged in vertical direction, lead- 
ing to a N x (NM) 2D network. 

(ii) We use a 2D tensor network of size M x M and di- 
mension D, arranged in the same plane as the initial 
2D network. One of the tensors has an open index 
corresponding to a physical system, while M 2 — 1 are 
auxiliary tensors. This leads to a (NM) x (NM) 2D 
tensor network depicted in Fig. [2;. 

(iii) We use a ID or 2D tensor network (see (i),(ii)), but 
arranged perpendicular to the initial 2D plane. This 
leads to a 3D tensor network as shown in Fig. [2ji- 

In each case, one may apply the method in an iterated 
fashion. For simplicity, we will consider only the networks 
at the first iteration. Similar as in the case of concatenated 
MPS, the states resulting in (i) can be interpreted as a 2D 
tensor network, but with increased (virtual) dimension x = 
D M in either horizontal or vertical direction. Similarly, in 
(ii) we obtain states corresponding to a 2D tensor network 
with a virtual dimension x = D M in horizontal and vertical 
direction. 

Note that already for very small M , the resulting states 
are useful, e.g., for a better approximation of ground states 
in 2D systems or to simulate time evolution in 2D. The 
advantage is that, while the underlying tensor structure 
is still two-dimensional - reflecting the geometry of a 2D 
system - one obtains with a relatively small overhead (a 
factor of M 2 ) an exponential increase of the (virtual) ten- 
sor dimension, Xcfi = D M . Given the fact that variational 
methods based on 2D tensor networks show a rather un- 
favorable scaling with the tensor dimension (0(D 12 ) for 
computational cost and 0(D 8 ) for memory [29j), one may 
use this approach to achieve virtual large tensor dimensions 



while keeping the dimension of the elementary tensors - and 
hence the computational cost and required memory - small. 
Although there is also an increase of the computational cost 
with the total number of tensors (scaling as 0(N 2 ) [29]); 
it is however much more favorable. Neglecting effects such 
as additional sweeps for optimization, the computational 
effort is increased by a factor of 0(M 2 ). For instance, if 
M = 3,D = 2, we obtain an additional computational over- 
head of roughly one order of magnitude due to larger num- 
ber of tensors, while the virtual tensor dimension is now 
given by 8. Using the initial TV x TV network with tensor 
dimension 8 would lead to a computational overhead factor 
of about 10 7 as compared to the D = 2 case. This differ- 
ence becomes even more drastic when considering larger D 
or M. 

The approach (iii) is the analogue of concatenated MPS 
we considered in the previous subsection. Similar as for 
concatenated MPS, one can describe all states resulting 
from a (post-selected) polynomially sized quantum compu- 
tation in this way if M = poly (AT). When considering a 3D 
tensor network as in (iii), one obtains single-qubit gates as 
well as nearest-neighbor gates acting on particles arranged 
on a rectangular 2D lattice. That is, the z-axis corresponds 
to the time axis, and the x — y plane corresponds to state 
of the TV x TV particles arranged on the 2D lattice after 
applying the (post) selected quantum circuit. Similar to 
concatenated MPS, also the interpretation in terms of a 
time evolution (of particles on a 2D array with nearest- 
neighbor couplings) is possible. Moreover, one may use the 
3D structure as ansatz states for a variational method to 
describe ground states or time-evolved states corresponding 
to some 2D systems. 



C. Tensor tree states with internal structure 

We consider now the example of a tree tensor network 
as shown in Fig. [3H- Tree tensor networks are quasi-one 
dimensional structures that can - similar to ID chains or 
MPS - be efficiently contracted in an exact way [HI]. I n our 
example, each of the tensors is of rank 3, and has dimen- 
sions di,d,2,ds, where d\ = d for tensors with open indices. 
We replace each of the tensors by a small tensor network, 
which we choose to be a triangle. That is, 

Di D 2 D 3 

— y ] y ] y ] B iia3a2 B i2a3ai B i3aia2 . (7) 

Ql — 1 Q2 — 1 CK3 — 1 

This process can now be iterated, i.e., each of the tensors 
A\ iakai is replaced by three tensors, say Cg^g 2 g 3 , m a tri- 
angular structure (see Fig. [3Jd)- There are two different 
types of tensors: External tensors - i.e., ones which are con- 
nected to outside initial tensors - of dimensions di, Dj, 
respectively, and internal tensors which have dimensions 
Di,Dj,D k . 

We consider a situation where Di < dj. In such a case, 
the internal structure of the initial tensor Ai t i 2 i 3 is deter- 
mined by elementary tensors C l g^g 2 g 3 , and in general this 
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FIG. 3: (Color online), (a) Graphical representation of a tree 
tensor network (TTN). The boxes correspond to tensors, where 
joint indices are summed over. Open indices correspond to phys- 
ical particles (red tensors), while auxiliary tensors are drawn in 
blue. The tensors are arranged in a tree-like structure, which 
guarantees that the contraction of the tensor network can be 
done efficiently, (b) Each of the tensors in the original tensor 
network is replaced by small triangular tensor structure in a 
concatenated fashion. 



restricts the values of Ai t i 2 i 3 . Notice that the entanglement 
features of the corresponding CTS, as measured by the en- 
tropy of entanglement, are determined by the dimensions 
of the tensors, and are in particular limited by the dimen- 
sion of the external links, i.e., d\,d2,d^. That is, in terms 
of entanglement, nothing can be gained by introducing the 
internal tensor structure. In order that the resulting tensor 
network state can carry the same amount of entanglement 
as the one described by the initial tree tensor network, one 
needs that the dimension of inner links at concatenation 
level k are larger than square root of the dimension of the 
links at concatenation level k + 1. In particular, D\ > y/di 
for k = 1, while for k = 2 tensor dimension (for the inner 
links) D% > y/D~i > d^ 4 are required. This can easily be 
seen by considering bipartitions of the system and by not- 
ing that the achievable Schmidt rank is determined by the 
dimension and the number of links between the two groups. 

The possible gain of such an internal tensor network 
structure is two-fold. First, the total number of parameters 
is reduced. While each initial tensor is described by cii^efa 
parameters, the resulting tensor network of depth k > 2 is 
described by [d\ + d2 + dz)D 2 +3 k ~ 2 D 3 parameters, where 
we assumed D k = D for all internal links. Second, the size 
of each of the tensors in the internal tensor network struc- 
ture is much smaller than the initial tensor. Many algo- 
rithms applied to the tensor network, e.g., the computation 
of normal forms of such tree tensor networks [31, 32, 33], 
or the optimization of tensors in a variational method [341 ] . 
scale with the dimension of the elementary tensors of the 
network. In spite of the usually polynomial scaling of these 
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algorithms, the computations quickly become intractable 
for increasing dk , so that a network containing tensors with 
small dimension are favorable in general. We have utilized 
this approach in [jjij], where numerical simulations using 
tree tensor networks are performed. 

We remark that the contraction of the resulting tensor 
network becomes more difficult as compared to the initial 
tree structure. This is due to the fact that the concatenated 
tensor network contains loops. To retain numerical accessi- 
bility, either approximate treatments have to be applied (as 
in contraction schemes introduced in the context of PEPS 
[HI) or the tree-like structure has to be kept, e.g., by lim- 
iting the tree- width of the concatenated tensor network (as 
in Ref. (3). 

IV. APPLICATIONS 

After having given some theoretical and analytical con- 
siderations for the possible advantages of CTS over other 
numerical methods for the description of states, we want 
to demonstrate applications of the CTS structure. The 
relevance of the CTS rests on two pillars. The first one is 
the ability of the (concatenated) tensor network to actually 
hold the relevant information about a state. The analytical 
considerations above indicate that this is the case for states 
based on circuits, time evolved states and others. The sec- 
ond pillar is the question if we can, once given a CTS, read 
out the contained information. Progress has been made 
with very similar networks in the context of PEPS. What 
we want to demonstrate in the following is the ability to 
find the potentially good description with numerically ac- 
cessible methods and see how good the approximation is. 
Moreover, this section has the aim to demonstrate the ap- 
plicability of the known contraction methods and describe 
some improvements thereof. 

A. The descriptive potential of CTS 

In this section, we want to demonstrate the descriptive 
potential of a CTS using a toy model. For reasons of com- 
parison, relevant states of the toy model were calculated ex- 
actly and these exact states were then approximated with 
both MPS and CTS. To not infiltrate the CTS descrip- 
tion with inaccuracies from an approximate read-out pro- 
cedure, we used an exact contraction algorithm for this 
network [H3|. 

In particular, we have tested the achievable accuracies 
when describing states resulting from time evolution in a 
spin chain, using the Hamiltonian 

H = ^i a) ^ +1) +B^\ (8) 

a a 

with B = 1 and a system size of N = 12 physical sites. 
The system is initialized in the product state 1+)®^ and 
evolved over a time T = 3.5, a point which is, in our units, 
close to the point where the fidelity of the CTS had a (pe- 
riodically recurring) minimum. Time evolution under this 



Hamiltonian shows the typical growth of entanglement in 
the state that makes MPS-based description hard. The 
optimal tensors in the CTS and also the MPS description 
were approximated by optimizing the overlap of the exactly 
calculated state and the tensor network state in a sweeping 
procedure. For each tensor, the overlap 

('tpex\'4>ex){'4>CTs\'ipCTs) , 

was calculated, leaving out one tensor to optimize. This 
tensor can then be found using linear algebra techni ques 
using the contraction result as a linear form. See Ref. [2a, 
123] and Appendix [Aj 

We compare the achievable accuracy when describing the 
state with MPS of varying dimension x> an d 2D CTS with 
varying numbers of rows of auxiliary tensors and tensor 
dimensions i.e., different M and Dk. These variations lead 
to the different number of parameters that the comparison 
is based on. Although a quadratic growth of the parameter 
count is expected for the MPS using this method, the plot 
shows a approximately linear growth. This is due to the 
fact that we did not count redundant parameters, which 
occur in the matrices close to the boundaries of the chain. 
We observe (see Fig. 2} that the description in terms of a 
CTS is more efficient, i.e., both a larger accuracy can be 
achieved when using the same number of parameters, and 
for a fixed number of parameters one can describe the time 
evolution accurately for longer times using CTS. 

Our tentative conclusion is that the additional structure 
leading to a reduction of the number of parameters and 
being imposed by the choice of CTS network reflects an in- 
ternal structure to be found in the time evolved state itself, 
comparable to the Suzuki- Trotter expanded time evolution 
operator that can be programmed into the CTS. This is 
supported by the interpretation of the rows in the network 
to be operators acting on an MPS (the very first row of 
tensors). It seems natural to assume that the rows con- 
tain a version of the time evolution operator of the system. 
However, these operators have not been programmed into 
the network this time, but found by the optimization algo- 
rithm alone. Further investigations are necessary, but the 
idea that an optimization algorithm together with a suit- 
able topological choice of network description yields a net- 
work of appropriate preparatory operators, reflecting deep 
structural properties of the described state, seems appeal- 
ing. 

B. Reading out information from CTS 

So far we have only considered the possibility to effi- 
ciently describe quantum states in terms of concatenated 
tensor networks, but not how to efficiently extract infor- 
mation from such a description or how to update it. Both 
the extraction of information and also updating procedures 
rely on the contraction of the tensor network, which is nor- 
mally used in a slightly modified version for this purpose. 

Given the fact that already a 2D tensor network is suf- 
ficient to describe all states resulting from a polynomially 
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number of parameters 

FIG. 4: (Color online). Comparison of the achievable accuracies 
when describing a time evolved state using MPS (blue, circles) 
and 2D CTS (red, squares). Simulated is the time evolution 
of a chain of N = 12 physical sites, initialized in the product 
state 1+)® and evolved over a time T = 3.5 using the Ising- 
Hamiltonian, Eq. ([8j, with B — 1. After this time T the chain 
is, in our units, close to the point where the fidelity of the CTS 
had a (periodically recurring) minimum. We compared MPS 
of varying dimension x> an d 2D CTS with varying numbers of 
rows of auxiliary tensors and tensor dimensions, leading to the 
different number of parameters. Shown in the plot is the value 
log(l — F), where F is the fidelity of the approximating state 
with the exact solution. Redundant parameters, as occuring in 
the boundary regions of the MPS, were not counted. We observe 
that the description in terms of a CTS is more efficient, i.e., a 
larger accuracy can be achieved when using the same number of 
parameters. 



sized quantum computation, one does not expect an effi- 
cient contraction of such a tensor network to be possible 
in general. In fact, it has been shown in [52j that con- 
tracting 2D tensor networks is computationally hard, the 
corresponding complexity class is #P. However, this does 
not mean that no efficient approximate methods can ex- 
ist which can successfully be applied in practice. In fact, 
in Ref. [lf| an approximate method to contract 2D ten- 
sor networks has been introduced and successfully applied 
in, e.g., the context of a ground state approximation for 
strongly correlated 2D systems [29]. This method will be 
described in the following, together with an investigation 
of two additional techniques: (i) A novel error correction 
scheme and (ii) An MPO compression scheme. 

For specific tensor network topologies (e.g., networks cor- 
responding to ID graphs, trees or networks with a bounded 
tree width [12]), an exact and efficient contraction and up- 
date of tensors is possible. If we are free to choose the 
contraction order of all indices, there exist Monte-Carlo 
based methods for the contraction, whose application will 
be shown in Appendix B. 



1. Approximate contraction of 2D tensor networks 

The approximate contraction of a 2D tensor network with 
open boundary conditions, as introduced in Ref. [26j], works 
as follows. The first (e.g., horizontal) line of tensors at the 
boundary can be interpreted as an MPS, where the lower 
indices are considered open. The second line can be viewed 
as a matrix product operator (MPO) acting on the first ma- 
trix product state. The resulting state (after contracting 
two lines) can again be described by an MPS, but of in- 
creased dimension. The aim is now to find (e.g., via a 
variational method) the optimal approximation of the re- 
sulting state by an MPS of fixed (low) dimension. This is, 
e.g., done by optimizing the individual tensors via solving a 
generalized eigenvalue problem (see Ref. f2(| or AppendixlAl 
below). The MPS found this way is now processed further, 
i.e., the MPO corresponding to the third line of tensors 
is applied, and one again aims at obtaining a proper ap- 
proximation of the resulting state by an MPS of fixed di- 
mension. The process is repeated until the second to last 
line of tensors is reached. The final step corresponds to 
calculating the overlap of the MPS resulting from above 
procedure (after processing all but the final line), and the 
MPS corresponding to the final line. All of these steps can 
be done efficiently. The evaluation of expectation values of 
(tensor product) observables works in a similar way. For 
details of the method, we refer the reader to j2§|. Notice 
that the same method can be used for 2D tensor networks 
where some of the tensors are auxiliary tensors (without 
open indices), as we consider in this paper. 

When using a concatenated MPS as described in 
Sec. IIII Al one may use the approximate method described 
above. However, especially when considering the descrip- 
tion of time evolution (Sec. IIII A 2|) or (post-selected) quan- 
tum circuits (Sec. IIII A ip . it is important to apply the 
method in a proper way, possibly utilizing symmetries of 
the state. In particular, in these cases the contraction 
should be done in the direction perpendicular to the or- 
der of the physical sites (left to right or right to left), 
rather than in lines parallel to the physical sites (up-down 
or down-up). A contraction in up-down direction would 
in these cases actually correspond to describing the state 
after each time step in terms of a fixed-sized MPS, and is 
actually equivalent to time evolution of an MPS as consid- 
ered, e.g., in Ref. [53]. When using a contraction in the 
perpendicular direction, such a limitation does not apply, 
see also [if}]. Numerical evidence suggests a significant in- 
crease in accuracy in this case. Moreover, the treatment 
of infinitely extended, translationally invariant states leads 
to the observation that a contraction over infinitely many 
columns of tensors perpendicular to the physical direction 
often results in a projection onto the eigenspace with the 
largest-magnitude eigenvalues of the MPO represented by 
the column. This makes it possible to employ additional 
exact numerical techniques, see, e.g., Ref. [161. 

If one is, like in the case of CTS, moreover able to choose 
the indices to contract freely, certain choices of tensors may 
allow for an efficient approximation via Monte-Carlo sam- 
pling techniques [38j, see Appendix B. There, the applica- 
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tion of Monte-Carlo methods to a 2D CTS will be demon- 
strated, using the inherent matrix product operator struc- 
ture of the CTS. We would also like to mention the possibil- 
ity to utilize String-bond state like tensor networks [39j in 
the context of CTS. Additionally, for certain choices of the 
tensors it is known that an exact and efficient contraction 
is possible [2l| . 

In the following we would like to suggest two improve- 
ments for the traditional contraction scheme. 



2. An error correction scheme 

We will now describe an error-correcting procedure for 
the contraction of 2D CTS, which is applicable also to the 
contraction of other rectangular grids including PEPS. 

We start with the traditional approximate contraction 
using the method described above, resulting in a number 
C, holding the contraction result. Following the line of 
argument from the sections above, we can interpret the 
number C as an approximation of the number 

C=(M 1 \M 2 ...M N _ 1 \M N ) 

where (Mi| is the MPS defined by the leftmost column 
of tensors, the operators Mj are the MPO defined by the 
columns in the middle and |Mjv) is the MPS defined by 
the rightmost column of tensors in the CTS. To remind 
the reader, a left to right contraction of the CTS involves 
the iteration of the following steps: (i) Start with i = 1 
and set (Mi| := (Mi|. (ii) Apply the MPO M t+1 to the 
intermediate MPS (Mi,...,i|. Both having a small bond- 
dimension, we obtain an MPS of large bond-dimension, 
(Mi,... ,i+i|. (iii) Reduce the bond-dimension of (Mi,.. ,; + i| 
to obtain another intermediate MPS (Mi,... ,j+i|, represent- 
ing (Mi,... ,j_|_i| as good as possible with this smaller bond- 
dimension, (iv) Increase i by one and continue with step 
(ii). The aim of the error correcting scheme is to estimate 
the error introduced by cutting off the bond-dimension of 
the intermediate matrix product states, and to correct the 
result C accordingly. 

More precisely, after the (i — l)th step of the standard 
left to right contraction, the CTS is approximated by 

C^Ci-i = (M 1 ,.., i \M i+1 ---M N - 1 \M N ) 

where (Mi,..,;| w (Mi|M 2 ---Mj. In the ith step we use 
the approximation (Mi,.. ,j_|_i| w (Mi,... ,j|Mj_|_i resulting in 

C « Ct = (Mi,.. , i+1 \M i+2 ■ ■ • Mjv-i|Mjv) 

The additional error of C in the ith approximation step 
is given by the value ej = Ci-i — Q, and the optimally 
corrected value of the contraction result is given by 

C = C + Y,ti- (9) 

i 

However, usually neither Cj_i nor Cj can be calcu- 
lated exactly since the exact MPS description of the state 



Mi + 2 • • • Mat_i|Mat) is too large to be computed. The cru- 
cial observation now is that the (N — i — 2)th step of a right 
to left contraction is a good approximation of this state 
with 

|M +2 ,.. iN ) m M i+2 ■ ■ ■ M N -x\M N ), 

which can be used to estimate the error 6j produced by the 
ith step of the left to right contraction 

6i = Ci_i - a (10) 
« (Mi 1 ..., i |M i+ i|M i+ 2,... 1 jv) 
-<Mi j ..., i+ i|M i+a ,..., lW ). 

This approximate value of is then used in Eq. [9l 

For an estimation of the achievable accuracy with this er- 
ror correction scheme, let the error of the overall left to right 
contraction be e. We note that also the error of the right to 
left contraction and its intermediate results (M,.. ,„| are of 
this size. Since the magnitude of the difference in Eq.fTOlis 
also of the order e, we are left with a residual absolute error 
of the order e 2 after the error correction. An application to 
toy models has confirmed our error estimation and yields a 
reduction of the error of about one order of e, or even bet- 
ter. For instance, the approximate contraction of the time 
evolved state in Sec. IIV Al with a cut-off bond-dimension 
D = 12 results in a value C with an error of 1.6(7)% with 
and 25(10)% without error correction, taking the mean of 
several approximations. 

While a similar reduction could in principle be achieved 
by using a bigger cut-off bond-dimension for the intermedi- 
ate results, the error correction scheme is favorable in most 
cases because of its better performance. As we can obtain 
all the required states (Mi,— , n | by caching one right to left 
contraction, we need merely twice the computation time for 
reducing the error by a factor of e. The overhead in mem- 
ory depends on the cut-off bond-dimension of the states 
(Mj,...,„|. Choosing this dimension equal to the dimension 
of the MPO Mi, the overhead is less than a factor of two, 
as we have to store N — 3 extra MPS which is less than the 
(N - 2)MPO + 2MPS of the CTS. 

We remark that the applicability of this error correction 
scheme is not restricted to CTS, but can in a similar way 
also be used e.g. in the context of the 2D PEPS approach. 

3. Compressibility of sequences of matrix product operators 

Additionally, the number of tensors in the CTS descrip- 
tion can be reduced significantly below the number needed 
in the canonical implementation of the Suzuki- Trotter pic- 
ture, as given in section fill A 2\ or for a generic network of 
(sparse) operators, like circuits. 

The reason is that it is not necessary to restrict each row 
to the description of a single Suzuki- Trotter (or generic op- 
erator circuit) time step only. Instead we can first put a 
good approximation for many of these rows, applied succes- 
sively, into one row, thus using the descriptive power of the 
CTS to the maximal extend. This is possible by calculating 
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and optimizing the overlap of one row of (variable) tensors 
with several concatenated rows of fixed tensors, in a way 
similar to maximizing the norm of a CTS when keeping ev- 
ery row but one fixed. We then concatenate these optimal 
rows, being fewer apparently, to reduce computational time 
in the read out process, whose computation time relies on 
the number of tensors involved. 

To get an idea of the potential of this ansatz let us con- 
sider time evolution. When performing time evolution by 
a Suzuki- Trotter expansion with MPO compression, there 
are two possible sources of error. The first kind of errors 
comes from the MPO approximation. This kind of error 
can be controlled, as we know the fidelity of the replace- 
ment step (the overlap of the rows to be compressed with 
the replacement row). When this fidelity is too small, we 
can reduce the number of rows to be compressed. The sec- 
ond kind of error comes from the Suzuki- Trotter expansion 
itself. This error can be made small by choosing a very 
small time step, so that the MPO corresponding to one 
row is close to unity. After compressing two rows to one, 
we are able to iterate the compression and compress two 
already compressed rows to one row, which now represents 
four time steps. This way, we obtain one operator cover- 
ing 2™ time steps with only n compression steps. Taking 
initially very small time steps thus does not result in a big 
performance hit, as the compression is very strong here, i.e. 
exponential, and compensates for it. The compression can 
of course only be applied as long as one row can in princi- 
ple hold the whole time step, but the results from the toy 
model in Sec. IIV Al make us optimistic that a row has the 
potential to hold comparatively big time steps. Numerical 
evidence concerning the compression fidelity for a variety 
of operators supports this view. In this case, the compres- 
sion not only improves the computation time but also the 
achievable precision by reducing the error introduced by 
Trotterization. 

The variational ansatz just shown leads to a picture of 
the CTS where a fixed number of tensors will be employed 
in an optimal way, as opposed to the direct programming of 
a set of analytically accessible operators into the network. 



4- Special cases - exact contraction 

Even though the problem of contracting an arbitrary 
2D tensor network is in general computationally hard (#P 
hard) , under certain conditions an efficient and exact con- 
traction of certain networks is possible. 

One such example are planar tensor networks, where 
each of the tensors fulfills a so-called match-gate (or free 
fermion) condition [2l[. It follows that if we restrict our- 
selves to CTS corresponding to planar structures, one can 
calculate the norm as well as expectation values of tensor 
product of observables efficiently for such states, as long 
as all tensors in the tensor network fulfill the match-gate 
condition. This implies that one may use such CTS, e.g., 
as variational ansatz states for ground states or time evolu- 
tion. In particular, we point out that the usage of auxiliary 
tensors as we propose for 2D CTS can be handled in exactly 



the same fashion. 

Another example are networks corresponding to trees or 
structures with a bounded or only logarithmically growing 
tree width. These also can be contracted efficiently and 
exactly. For instance, the contraction of a subcubic ten- 
sor tree (i.e., a tensor network where each of the tensors 
in the tree is connected with three or less neighbors) has 
a computational effort scaling as 0(D 3 ). For a variational 
method for the search for ground states and the maximiza- 
tion of overlaps with CTS enhanced tree tensor networks, 
see Appendix A. 

5. The advantage of CTS in the operator picture 

Being able to program the Suzuki- Trotter expanded time 
evolution operator or other generic quantum circuits di- 
rectly into the state description offers access to alterna- 
tive advantageous numerical approaches. Time evolution 
methods usually rely on maximizing a significant number 
of overlaps of the kind 

\(A+sM5t)\ipt)\ 2 

(lpt+St\ll>t+8t) {iptHt) 

where r(St) is the time evolution operator for a time St, \ipt) 
is a known tensor network state and the tensors describing 
the state \^ t +8t) have to be found. In the context of cir- 
cuits, an application of a set of gates can be regarded as a 
time step like above. Starting from this expression, to com- 
pute the time evolved state after a time T, one would have 
to compute single time steps repeatedly, and each time one 
would have to perform network contractions to determine 
the optimal tensor entries. (Algorithms of this kind are 
found to converge to a reasonable approximation of the best 
tensor network description of the desired time evolved state 
\i/>t))- Depending on the implementation, finding the opti- 
mal tensors can consist of many sub-steps, e.g., a sweeping 
procedure approximating single tensors while leaving the 
remaining tensors fixed, each sub-step requiring another 
contraction. 

The CTS description of states is very efficient in this 
regard if, as given above, the time evolution operator 
r (^) = I1 T (^) ls programmed into the structure and de- 
scription of the state itself. Using CTS, we are able to 
circumvent the many contractions and possible sweeping 
steps by a single optimized contraction. In contrast to a 
traditional MPS time evolution, for example, it is possible 
to use error correction and alternative contraction order 
(e.g., left-to-right instead of top-to-bottom). Moreover, we 
are able to employ MPO compression, which is, however, 
also applicable to the traditional time evolution method, 
but there not in a direction-optimized fashion. 



V. SUMMARY AND OUTLOOK 

We have introduced concatenated tensor network states 
(CTS), a class of states that is obtained by decomposing 
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the high-rank tensor describing the coefficient of a multi- 
particle states into a tensor network in an iterative fashion. 
The basic idea is to impose additional structure to each of 
the tensors appearing in a tensor network description of a 
given state. We have demonstrated this approach for ID 
tensor networks, where in a first step a description in terms 
of a matrix product state is obtained. Each of the matrices 
(tensors) is then further decomposed into a ID tensor net- 
work (matrix product operator), yielding a 2D tensor net- 
work with many auxiliary tensors in the next step. Similar 
methods can be applied to 2D systems, yielding 2D or 3D 
PEPS with auxiliary tensors, or to tree tensor networks. 

We have demonstrated that with such CTS, one can de- 
scribe multi-particle quantum systems with rich entangle- 
ment features in an efficient way. In particular, states aris- 
ing from time evolution or generated by polynomial (post- 
selected) quantum circuits can be described, and an inter- 
pretation in terms of (post selected) measurement-based 
quantum computation can be given. The states can - in 
contrast to matrix product states or projected entangled 
pair states - contain a large amount of block-wise entan- 
glement (up to a volume law) and long-ranged correlations, 
while their description remains efficient. In particular, a 
subclass of matrix product states and projected entangled 
pair states with high effective bond dimension can be de- 
scribed. 

We have demonstrated that it is possible to describe 
states arising from time evolution of a ID quantum sys- 
tem with help of such a 2D CTS more efficiently than with 
a matrix product state. We have discussed the description 
arising from a Trotter decomposition of the evolution oper- 
ator, as well as direct optimization of (auxiliary) tensors in 
the 2D tensor network of given size and dimension. In this 
context, we have also applied a method to compress matrix 
product operators to obtain a more efficient description of 
the time-evolved state. 

We have also discussed and improved methods to read 
out information from 2D tensor networks. The applica- 
bility of approximate contraction methods, possibly with 
different direction of contraction (left to right), has been 
discussed and improved using an error correction scheme. 
Both the effectiveness of the CTS description in the context 
of time evolution of one-dimensional systems as well as the 
impact of our suggested enhancements to the traditional 
read-out methods were demonstrated using numerical re- 
sults for a toy model. Also the applicability of Monte-Carlo 
methods for the contraction was demonstrated. 

The results indicate that the new class of states is use- 
ful in the context of describing and simulating time evolu- 
tion of ID quantum systems, but might also be used for 
the simulation of ground states of 2D quantum systems. 
The different interpretations in terms of trotter decompo- 
sition, (post selected) quantum networks or (post selected) 
measurement-based quantum computation we provide may 
also inspire a new point of view to tensor network states 
and encourage further development. 
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APPENDIX A: VARIATIONAL OPTIMIZATION 
OF CTS-ENHANCED TREE TENSOR NETWORKS 

Important applications of quantum mechanical simula- 
tions are the search for ground states and the computation 
of the time evolution of states governed by a given Hamilto- 
nian, usually employing variational methods. On the math- 
ematical level, an essential element of the variational proce- 
dures in tensor networks is the linear dependence between 
the network and each of its tensors. Contracting the CTS 
(or an amplitude or expectation value involving a CTS) 
leaving out one of the tensors provides us with a simple 
linear or quadratic form which is suitable for investigation. 
The maximization of an overlap or minimization of an en- 
ergy is thus reduced to the analysis of such a form and can 
be performed using linear algebra. Naturally a possibly 
exact and efficient contraction method is desired. 

More precisely, for instance, finding the ground state of 
a Hamiltonian H means finding the state that solves 



H\i/>) 



mm. 



We can write a tensor network state as 



\TNS) 



^ t ^si...s„ai. 



l Rs n+1 ...s N a 1 ...a n I s ) 



(Al) 



where T is the tensor under consideration and R is the 
remainder of the tensor network, already contracted up to 
the indices that connect T and R. The mentioned linear 
dependence on T is exploited by (virtually) replacing the 
tensor T by tensors D (s, a) which have the entries 



1 s = s and a = a 
else 



(A2) 



With help of the tensors D we generate states 

IV' (s, a)> ■=^D{s 1 k) Si Snai _ an R Sn+l ... SNai ... an |s) 

s.a 

where (s, a) is a combined index. With these states we in 
turn generate matrices 



as well as 



%,&),(§',&') := <^(s,a)|iJ|V'(s / ,a')) 



N(s,a),(s',a') := (i>(s,a)\ip (s',5.')). 
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Finding the entries of the tensor T is now reduced to a 
generalized eigenvalue problem, 



t* -E-t 



mm, 



t* ■ N ■ t 

where the tensor t with the minimum generalized eigenvalue 
E ■ t = \N ■ t 

is the solution of the local minimization problem, i.e., the 
minimization problem with respect to T when the other 
tensors are fixed. The exact technical implementation of 
this idea is, of course, subject to optimization and will not 
be done by the mentioned contraction over dummy tensors. 

We are facing two numerical problems. One is the con- 
traction of the tensor network. The second problem is 
finding the generalized eigenvalues of the matrices given 
in the section above. The difficulty of the contraction of 
a network increases polynomially with the index rank and 
moreover depends strongly on the topological structure of 
the network. Tree networks can be contracted with efficient 
algorithms [3l[ , and the dimension of the tensors is the pa- 
rameter which governs the efficiency of the contraction in 
this case. Contracting networks with loops on the other 
hand is in general intractable if performed exactly. 

Obviously, the nested (e.g, triangular) tensor structure 
does not simplify the contraction of the network, but re- 
mains feasible if the tree-width of the tensor network is 
small. A tree whose tensors are replaced by small loops is 
such a tree-like structure. In the corresponding efficiency 
considerations, the role of the dimension connecting the 
tensors is replaced by the dimension connecting the loops 
among each other. Now, there are values of the dimen- 
sion (entering polynomially into the computational effort) 
where a contraction of the network is still possible, but a 
solution of the generalized eigenvalue problem is not - the 
reason being the size of the corresponding matrices, whose 
size is scaling like D 3 x D 3 if D is the index rank. 

The nested tensor loops address this problem by "shield- 
ing" the large outgoing index. More precisely, let us define 
a tensor network state like in Eq. IA11 For the sake of sim- 
plicity we consider a tensor A not connected to any physical 
sites and the network to be a subcubic tree, 



\TNS) :=2> Q1 



l + 1 ...s N a 1 a 2 a 3 



We are now able to rewrite the state \TNS) by replacing 
A by a loop as shown in Fig. [3] 



A ai a 2 a 3 '■— y ] ^aiQ(3-^a 2 a7^. 

a/37 



3 

03/37- 



where now correspondingly 



\TNS loop 



EE 

s,a a/37 



rjl rj2 r>3 d 

ai a/3 JD a 2 a-f JD a 3 f3'y rl s„+i...SNai 0203 



This helps to reduce the size of the matrix of the cor- 
responding eigenvalue problem because a) the tensors 
B^.B^.B 3 can be locally optimized individually, and b) 
the indices a,/3, 7 can have smaller dimension, while the 
loop structure of the tensor T replacement network retains 
the entanglement properties which are so important for the 
power of the description. It is possible to choose low but 
sufficiently high index rank for the internal indices such 
that the entanglement being carried by the external indices 
(that connect the loops among each other) is not reduced. 

In detail, finding the optimal values of the loop tensors B l 
can be performed as follows. First, the network represented 
by R has to be contracted. Once this tensor is found, it 
is kept fixed for the optimization of the tensors B % . We 
then repeat the optimization steps for the loop tensors as 
described in the section above, using the state 



? 3 

03/37 



\TNS loop ) := E E dl fa' fi ' „ B l. ai B 3 a 

— ^7 \ ' aia/3 

X Rs rl+1 ...s N a 1 a 2 a 3 | s ) > 



a/37 



with a tensor d 1 like in Eq. IA2I Similarly, we proceed for 
the tensors B 2 ,B 3 . In these steps we can make use of the 
fact that several (more than one) sweeps through the loop 
tensors will give a better convergence, while the computa- 
tional overhead for this is small, because the huge remain- 
der of the network - represented by the tensor R - stays 
constant and does not need to be contracted again. If the 
dimension of the internal indices is large enough, several 
sweeps through the loop will converge to a network 



E B a ia fS B . 

a/37 



2 R 3 



with a tensor A whose values are the same as in the case 
without the replacement network. In some cases the origi- 
nal problem of finding A would not have been feasible, but 
even if so, the sweeping through the loop gives an advantage 
in computation time. 



APPENDIX B: MONTE-CARLO SAMPLING OF 

CTS 

In some instances, it is possible to contract the concate- 
nated tensor network approximately with a Monte-Carlo 
based approach. For this, let us quickly recall how the 
Monte-Carlo method works. The easiest and most basic 
Monte-Carlo (MC) technique is the Metropolis algorithm 
(37| . Like all MC methods, it is used to estimate inte- 
grals (or sums) over high dimensional integration spaces. 
In these spaces, naive approaches like Riemann-integration 
require a huge number of sampling points for a certain re- 
quired accuracy, whereas usually the MC methods show a 
much quicker convergence to the exact value. 

The basic idea is that we can select a sample of points in 
the integration space such that 



f (x) n (x) dx re AT 



E 



At;, 
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where Av is a unit volume in V, /i is a well-behaved mea- 
sure on V, and Z — J v /i (x) da;. Naturally, the selection 
rule for the set of samples, {xi}, is the key and has a foun- 
dation in statistical mechanics. Assume that / is a prop- 
erty of an ergodic physical system with density (probability 
density to be found at that point) n in configuration space. 
The system being ergodic, we know that the time average 
of the property / equals the average of / over configuration 
space with weight /i, 



(f)t = / f(x)fi(x)dx. 
Jv 

We obtain the set of samples {x{\ by simulating the be- 
havior of the system in time and recording the position 
x{U) = Xi at discretely (and equally) spaced points {U} in 
time. Let now P (x — > a;') be the probability of the system 
to go, during one discrete time step of a random walk, from 
point x to point x' . A set {x^ of a random walk derived 
with such a rule is called a Markov chain, with the essential 
property being that Xi is only dependent on (and not 
Xi-2 etc.). It is known that the so called detailed balance 
condition for the probability P, 

H [x) P(x^ x') = fi (a/) P (x' -> x) , 

is a sufficient criterion to ensure that a random walk of 
the system, ruled by the transfer probability P, yields a 
time average approaching the value Z^ 1 J v f (x) [i (x) dx 
for t — > oo. One can impose this transfer probability by 
the following rule: 

1. Being at point Xi, choose randomly a position £. 

2. Calculate the value A (x+ — > £) = min (l, j^§7) J ■ 

3. Randomize a number in the interval a G [0, 1]. 

4. If a < A(xi — > £), then Xi+i = £, otherwise Xi+i = 

X i . 

This is the (basic) Metropolis algorithm (37| . The prob- 
ability of going from x to a;' under this algorithm obeys the 
detailed balance condition and hence yields a sample that 
is representative for the measure fi. We note that with this 
rule we can generate arbitrarily large sets of positions in 
time without the need to store the set {x{\ itself. Fur- 
thermore it is not necessary (for the evaluation of the time 
evolution) to know the value Z = J v fx (x) dx, which cancels 
in the calculation of A; a fact that makes it possible to work 
with relative probabilities and unnormalized measures. 

We now want to show that the contraction of the con- 
catenated tensor networks can be implemented via a MC 
algorithm. To demonstrate the principle we give the formu- 
las to contract a toroidal network of N x M tensors of rank 
4, although the formalism is easily adapted to non-toroidal 
networks and higher dimensions. Consequently, we want to 
calculate 



E 

indices s„,s d ,s,,s r 



IK 



where u,d,l,r mean "up, down, left, right" respectively, 
the indices s Ui d,i, r depend on the position within the 
network, and s d = s u (i + l,j), s r = s t (i,j + 1), 
s u (1, i) = Sd (N,j) and si (i, 1) = s r (i,N), 

The basic principle is to perform the contraction over 
the indices s Ut d,i,r in a hierarchical order: We first contract 
over the indices in each row. This is formally the trace 
over a product of matrices, the matrices being the tensors 
of rank four, where the indices connecting in the vertical 
direction are kept fixed. In the next step, we contract over 
the indices that connect the rows. Following this idea, in 
the case of an n-dimensional network, the hierarchy has n 
levels - indices of increasing level thereby connecting slices 
of increasing dimensionality. For the problem at hand, we 
write 



j)»l(i,j)si(ij + l) 



Su,Sl 



En^.(i>.(i+i)' ( bi ) 



where s u (i) = (s u (i, 1) , s u (i, 2) , ...), 



B« (<>«(*+!) 



Tr 



l[T^[s u (i,j),s u (i + lj)] 



and where T*'- 7 



are matrices with elements 



(TU[a,b]) t 



Su(i,j)sd(i,j)si(i,j)s r 



What we see is that, while keeping the up and down in- 
dices fixed in each row, we obtain a set of matrices for each 
row, and carry out the summation of the left and right in- 
dices with matrix products and a trace. As seen in Eq. IB1| 
this leaves us with another set of matrices, R l s 

which are treated like the matrices T 4J [•,••]: they are mul- 
tiplied and the product is traced over. This means that the 
contraction over the whole set of indices is the trace over 
a matrix product for matrices whose elements are again 
traces over matrix products. This recursion repeats itself 
and adds another "generation" of matrix products for every 
dimension of the tensor network to be contracted. If the 
tensor network was no torus, the boundaries would have 
been taken by tensors of lower rank. This would then re- 
place the trace by the contraction with the tensor at the 
boundary; the situation is similar to matrix product states 
with periodic versus open boundary conditions. 

As we see, formally, we use matrix products and traces 
over all levels of matrices in the hierarchy, but for the sec- 
ond generation of matrix products and higher ones this 
cannot be carried out explicitly anymore. While the di- 
mensions of the T IJ [-, ••] -matrices depend on our arbitrary 
choice of the index rank of the tensors T only - making 
the first generation of matrix products feasible - the di- 
mensions of the second generation of matrices (here: the 
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row matrices R) is already much too large: The number of 
columns enters exponentially. This is the point where the 
Monte-Carlo sampling is used. We sample the (second and 
higher) generation of matrix products. 

To improve readability, we leave out the index u from 
now on and write the matrix elements in the following way 
R l s , i+1 s = (si \ R l \si + i) (as usual in quantum mechan- 
ics). There are in principle two possible ways how to per- 
form the MC-approach in detail. 

In the first, the contraction takes the form 



E 

Si,...,S]V + l 



S 1 |i? 1 | S2 )(s 2 |i? 2 ...ir|Sn+l) 



x (s„ +1 | R n+1 |s n+2 ) (s n+2 | R n+2 ...R N \s N+1 ) 

= E L ( s l, ■■■, s N+l) R(S1, ■■■,S N+1 ) 



si....,sjv + i 

= E |i(si,---,s«+i 

si,...,sjv + i 

= E 

Bx,...,S N+1 



,2 R (s n +l, Sjy+l) 

L (si, s, i+ i)* 
)f( Sl ,...,s N+1 ). (B2) 



where 

L( Sl ,...,s n+1 ) := (s 1 \R 1 \s 2 )(s 2 \R 2 ...R n \s n+1 ) 

R ( s n+l, s N+l) '■= ( s n+l| R" +1 |s„+ 2 ) 

D JV i„ 

1+1/ > 



x (s n+2 \R n+2 ...R N \s N+1 ) 



L and \i are actually independent of s/., k > n + 2 and R is 
independent of a*, k < n. This way, we have derived a for- 
mal expression that is compatible with an MC-approach: \x 
is a positive measure (obviously, fj, > 0), / is a function to 
be integrated over and weighed by /i. The integration (sum- 
mation) space is the set of all vectors S = (si, ...,sat+i). 

The second, alternative, implementation of the MC- 
approach is similar, but divides the integration space into 
smaller subsections. We use the same ansatz, but this time, 
we do not look at the functions L and R, but at the func- 
tions I and r 



^( s li s ri+l) ; — 2_^i ^ ( Sl ' '"' S «+l) 

S2,...,Sti 

r(s n+ i,sjv+i) := i2(s n+ i,...,sjv + i) 



instead. Then 

En^oosfi+i) 



E,,/ 2 r(s n+1 ,s N+1 ) 

/ (Si,S„+i) 



Sl,Sn+l,SjV + l 



where the sum is performed over only three (instead of 
N+l) vectors. The functions I and r, which are sums them- 
selves, can be calculated, again, with an MC-approach. The 



expression 

i(si,S B +l)= Yl (S1\R 1 \S2){S2\R 2 -R n \s n+1 ) 

S2,-..,S n 

(and likewise r) are very similar formally to the original ex- 
pression for the whole sum in Eq. lB2l and can be treated in 
an analogous fashion: another subdivision of the summa- 
tion space {s 2 , s„} into two parts will be initiated. This 
recursion will have a depth of order log (N); hence, assum- 
ing we need for a summation over three vectors O (m) MC- 
samples, this method will lead to a computational effort of 

O ((2m) logAr ) = O (2 l °z N N y °z m ). The decision to choose 

way 1 or way 2 should be based on tests concerning the 
achieved convergence speed and depend on the problem at 
hand. 

There is one more bookkeeping-issue to consider in both 
approaches. Let & [•] denote the limes of the Monte-Carlo 
sampling over a Markov chain {Si}, approaching infinite 
length, that is generated by the measure [i. Then 



E^ 1 



6[/] 



all S 



E^( S )/( S ) 

all s 



that is, the Monte-Carlo summation yields a result as if 
the measure was normalized. Hence we need an estimate 
for the partition sum Z — X^all s A* (S)- It can be obtained 
through the relationship 



E 

all s 



A1 (S) / i- 1 (S)=Z- 1 |{S}| 



hence Z = |{S}| (S [m -1 ]) • This doubles the computa- 
tional effort. 

The applicability of the MC method is not universal 
though, the limits being set by the measure \i and the func- 
tion /. As mentioned in the description of the MC method, 
the system needs to be ergodic. This is not the case if the 
measure fi has regions whose borders cannot be crossed 
during a random walk. This phenomenon is usually due to 
the measure being zero (or very small compared to other 
regions) in some regions of the configuration space. In this 
case, the random walker decides too often to stay where it 
is and traverses the configuration space either too slowly or 
not at all, obviously not yielding an interesting sample col- 
lection. A property that can cause non-ergodicity is hence 
sparseness of the tensors in the concatenated tensor net- 
work. A remedy for this problem is a local change of basis, 
and hence readily available. 

Another, more serious, problem is non-convergence of the 
MC sampling due to the so called sign-problem. This means 
that the relative error of the sampling method stays big 
over an exponentially large sampling set, because, roughly 
speaking, in the summation process two very big sub- 
summands (with small relative error) cancel each other. 
This problem does not occur if all (row-)tensor elements 
are > 0. 

Also other tensor structures are known to be accessible 
in similar ways. For instance, so called string-bond states 
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[39| can be efficiently contracted via Monte-Carlo sampling. 
While in the initial work only string-bond states where all 
tensors are associated with physical particles have been 
considered, one can easily adopt the construction of con- 
catenated tensor network states to string-bond states. In 
particular, can use exactly the same tensor structures as for 
string-bond states, where however only some of the tensors 



correspond to physical particles, while the other tensors are 
auxiliary. 

We have performed numerical tests on this method and 
found that for certain cases of tensor networks, e.g., with 
positive tensor entries, the algorithm converges quickly to 
a value of the contraction result with a quickly decreasing 
relative error. 
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